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We derive the boundary conditions for M 0 S 2 and similar transition-metal dichalcogenide hon¬ 
eycomb (2H polytype) monolayers with the same type of fc-p Hamiltonian within the continuum 
model around the K points. In an effective 2-band description, the electron-hole symmetry break¬ 
ing quadratic terms are also taken into account. We model the effect of the edges with a linear 
edge constraint method that has been applied previously to graphene. Focusing mainly on zigzag 
edges, we find that different reconstruction geometries with different edge-atoms can generally be 
described with one scalar parameter varying between 0 and 27r. We analyze the edge states and their 
dispersion relation in M 0 S 2 in particular, and we find good agreement with the results of previous 
density functional theory calculations for various edge types. 

PACS numbers: 73.20.At,73.21.Hb,73.22.-f,73.22.Dj,73.61.Le,73.63.Bd 


I. INTRODUCTION 

Semiconducting transition metal dichalcogenides 
(TMDCs) have attracted much attention lately. This 
applies especially to their monolayers, which have a 
direct band gap (M0S2, MoSe2, MoTe2, WS2 and WSe2) 
as opposed to their three-dimensional crystals with an 
indirect band gap. TMDCs are interesting for a number 
of applications ranging from flexible and transparent 
field-effect transistorJ^H to logicaP and optoelectronic 
devices These applications necessitate the knowl¬ 
edge of the electronic structure of finite flakes and 
nanoribbons, which in turn also relies on the accurate 
description of the edgesStudying the edge states is 
also very important when charge carriers need to be 
injected into a layer of TMDC. The edges of TMDCs 
have alrea dy b een studied experimentall^y^^^ and 
theoreticall}Ml^M^ by several groups. Recently, using 
STM spectroscopy, experimental evidence of metallic 
edge states has been found in single layer MoS2f^ 

Ah initio calculations, like density functional theory 
(DFT) combined with tight-binding (TB) modelJ^have 
allowed many detailed studies of TMDCs. Nevertheless, 
the description of mesoscopic edge-terminated systems 
still represents an outstanding challenge for theory. In 
such cases, a continuum model can offer an efficient way 
to describe the system. Here, we study the boundary 
conditions (BCs) and the low-energy in-gap edge states 
in TMDCs in general, and in M0S2 in particular. In¬ 
spired b y prev ious works on the analogous problem in 
graphenewe derive the general boundary conditions 
for the zigzag and armchair types of edges in monolayer 
TMDCs, see Fig .[2 The BCs are given in the form of two 
linear matrix equations, one for the wave function and 
one for its derivative normal to the edge. These equa¬ 
tions can be expressed using an initially unknown matrix 
M, in which the number of free parameters is reduced 
using the symmetries of the system and the edge. Due 
to the continuum nature of our method, not all the pa¬ 
rameters in M can be determined exclusively by means 


of symmetry considerations. Some atomistic details of 
the boundary elude this procedure, and therefore, some 
parameters need to be fitted to results obtained from 
other methods with atomistic resolution, such as DFT or 
density functional based tight-binding (DFTB) methods. 
We will demonstrate this fitting procedure through the 
example of zigzag edges in M0S2. 

To describe a finite system in M0S2, we start with the 
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FIG. 1. The honeycomb lattice of transition-metal dichalco¬ 
genide monolayers. Rectangles with dashed lines show the 
unit cells of the armchair and zigzag nanoribbons. The width 
of the nanoribbons can be parametrized by the number Nd 
of the A-B dimers in the corresponding honeycomb unit cell 
without regarding any adatoms attached to the possibly re¬ 
constructed edge. In this sense, the width of the framed arm¬ 
chair ribbon would be 14, and that of the framed zigzag ribbon 
would be 7 in this figure. The actual width of an armchair 
ribbon is ^Nda, and that of a zigzag ribbon is ^Nda, where 
the lattice constant in M 0 S 2 is a = |ai| = \a 2 \ = 3.1565A.^ 
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bulk Hamiltonian for one spin state, 
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_ ( 1 ) 

In this form, the effective k p HamiltoniarP^ is cast into 
the valley isotropic represent at ion,where the blocks of 
H correspond to the K and K' valleys, respectively. The 
basis functions, consisting of mainly the d-orbitals of the 
Mo atoms, are the wave functions of the valence band 
and the conduction band, as in Ref. m- The basis of 
the 4 -spinor thus reads 

valley, the valence band maximum is £y\ the conduc¬ 
tion band minimum is £c ^, and we also take into account 
the quadratic terms with and which break the 
electron-hole symmetry. We will refer to the nonequiva¬ 
lence of the two blocks of H, which is due to the spin-orbit 
interaction (SOI), as valley asymmetry. We assume that 
7 is the same in both valleys, and since we focus on the 
vicinity of the K points, we neglect the trigonal warping 
and cubic terms.^The distance from the K points is de¬ 
noted by k, and k± = kx ^ iky. The eigenvalues of the 
block of H for the R" valley are 


£± — - (^v + 5c + k‘^{a + /3)) 


± \lk‘^-f‘^ + - {k‘^{a -p)- E^Y 


( 2 ) 


while £'±{k') in the K' valley can be obtained with the 
respective parameters. The energy gap between the 
conduction- and valence-band edges is . 

The corresponding eigenstates are 


5 T E„l‘ 2 i — c^k^ 


0 , 0 , 1 , 


^k- 

-E'J 2 ^ 


, 0,0 


p'k'' 


^k'_ 


in the K, and 


in the K' valley, ( 3 ) 


where T denotes transposition. At an edge, where only 
the parallel component of the wave vector k\\ remains 
a good quantum number, the degeneracy is eightfold: 
twofold for the two valleys, twofold for the two solutions 
of the wave number at a given 5 = = £^_!}{k2^) 

energy, and twofold for the ± sign of the wave vector’s 
normal component k_\_^ see Fig. For now, we neglect 
the electron spin, and we assume that the armchair and 
zigzag edges mix these eight states to satisfy the bound¬ 
ary conditions. 

The Hamiltonian H is not specific to M0S2, as it 
can also describe similar transition-metal dichalcogenide 
monolayers with honeycomb lattices, like MoSe2, MoTe2, 
WS2 and WSezPxhe M matrix method that we present 
in this paper is also as general as the Hamiltonian itself. 
Nevertheless, in the last section, we will demonstrate the 
use of our method through the example of MoS2p^ 



FIG. 2. Schematic representation of the eight degenerate 
states in /c-space that are mixed by the armchair and zigzag 
edges. These in-gap states denoted by dots are characterized 
by the same spin, the same energy e and the same quantum 
number /cy, where k\\ is the component of k parallel to the 
edge. The dispersion curves are labeled with the correspond¬ 
ing functional forms, which are the diagonal elements of the 
Hamiltonian H. For solutions localized at the edge, the nor¬ 
mal component k± is imaginary, which results in the opposite 
curvature of the valence and conduction bands. States from 
these two bands with the same wave vector compose the cor¬ 
responding in-gap edge state. From the left to the right, the 
normal components of the wave vectors are labeled as —ik±i, 
—ik±2, ik±2, ik±i, —ik'^i, —ik'^2, ik'^2 and ik'^i. 


II. THE M MATRIX METHOD 

Compared to the calculation for graphenej ^^ l ^^ * an 
additional boundary condition is needed here for the 
first derivative of the wave function T', because of the 
quadratic terms in H. As we will see, the second con¬ 
dition can also be expressed by the same, constant M 
matrix. Inspired by Ref. we introduce hard wall 

boundary conditions for T and T' with the hermitian 
matrix, Mq: 


[u + Mo 5 {n-r)]^{r)=E^{r), ( 4 ) 

where 6 {n'r) is the Dirac delta function modeling the 
confinement potential at the boundary, which crosses the 
origin here for simplicity, and the unit vector n is normal 
to the boundary and pointing outwards. To obtain infor¬ 
mation about the unknown matrix Mq, we first consider 
the infinitesimal integral of Eq. @ through the bound¬ 
ary: 


[h^ + + Mo5{n-r)]^{r) n-dr = 0 , (5) 

where the Hamiltonian is split into parts with zeroth or¬ 
der, linear and quadratic terms in /c as H = H^;^ -f- + 

Hq,,/ 3 with indices referring to the parameters they de¬ 
pend on. The integrated zeroth order part and the right 
hand side of vanish for any finite wave function. Af¬ 
ter carrying out the integral in the given limit, we arrive 
at the following equation for T = '^{r)\n-r=o- fhe 
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boundary: 

7J^5' + A^' + ^Mo^' = 0, (6) 

where = ^^(^)|n r=o- denotes the derivative nor¬ 
mal to the boundary, the diagonal matrix A = 
diag(a, /3', Of'), and = tq 0 exp (i|n-cr). Here and 
throughout the paper s, r and a are denoting vectors 
of the Pauli matrices acting in the spin space, valley 
space and the space of the basis functions ipc}^ with 
the zeroth component being the two-dimensional iden¬ 
tity matrix. Note that the operator is unitary and 
antihermitian, and so = —1. Moreover, A combined 
with the derivative operator is antihermitian. Therefore 
if we calculate the expectation value of these terms by 
multiplying them by T”*" from the left, the antihermitian 
operators yield purely imaginary numbers, but the her- 
mitian matrix Mq yields a real number, meaning that 
they have to be zero separately: 


the case. Therefore, the number of the linearly indepen¬ 
dent solutions for T must be less than four. One could 
argue that T' should be treated as four additional in¬ 
dependent variables adding four more dimensions to the 
Hilbert space, opening up the possibility to solutions that 
are linearly independent in the subspace of T', but not 
in the subspace of T. This option can be dismissed if we 
regard Eq. which fully determines T' for a given T, 
because A is invertible. This means that the subspace 
of the solutions in T will determine the corresponding 
subspace in T' with the same dimensionality. 

We now move to the basis spanned by the solutions to 
T, and we denote the size of this subspace by d, which 
is 1 < d < 3. The first d coordinates of T can thus be 
chosen arbitrarily and the last 4 — d have to be zero. For 
Eq. (12), this means that the top left d x d block of Ty 
must be 0 in this basis. At the same time, the last 4 — d 
columns of Mq can be chosen freely. We then define the 
matrix Mi in the following way: 


74/0^4/ + = 0 and 

(7) 

( 0 

Afi — 1 

-yt' 

= 0. 

(8) 

IVll — 1 

V-^ 

0 


(14) 


A second consideration concerns the total current, 
which can be determined as 


J = |[H,n-r] = 
/ —2ian-'V 
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ye 

—2i/3n-V 
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—2i/3'n-V 
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(9) 

0 \ 
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ye 


-*x 


—2m'n-V J 


where T is a (4—d) x d-sized matrix containing the matrix 
elements of the bottom left block of (row indices from 
d + 1 to 4 and column indices from 1 to d). The top 
left d X d-sized block contains zeros, and the same holds 
for the bottom right (4 — d) x (4 — d)-sized block. Since 
Ty is unitary and antihermitian, it is easy to prove that 
it consists of the same blocks as Mi, except for the sign 
change of the bottom left block. It follows from here that 


measured in the direction of n = (cosy, siny), which is 
normal to the boundary, and V stands for the gradient 
operator. We split J in diagonal and off-diagonal parts: 


7 

J = -\- Ja 6 = —illTo (8) exp 

h 

( 10 ) 

where Ty = It is a natural requirement as a 

boundary condition that the expectation value of the nor¬ 
mal current at the boundary needs to vanish: 


(J) = -i|(j^)-i^(^(n-V))=0, (11) 

which can be translated into yT’t'TyT + 2T’t'AT' = 0 at 

the boundary. Comparing this to Eq. ([^, we see imme¬ 
diately that 

= 0 and (12) 

TUT' = 0, (13) 


.TT 

^2 


\ 2 

ri'Crj — ij^A{n • V), 


meaning that the expectation values of the diagonal and 
off-diagonal currents both vanish: {Ja,p) = (d^) = 0- 
We assume now that the boundary problem has at least 
one solution. From Eq. (12), it follows that we could 
only have four linearly independent solutions in the four¬ 
dimensional Hilbert space of T, if Ty was 0, which is not 


TyT + MiT = 0 


(15) 


for every solution of the boundary problem, and Mi is 
a hermitian matrix. Eq. (15) multiplied by y and sub¬ 
tracted from Eq. leads us to 


A^' + M 2 4/ = 0, 


(16) 


where we introduced M 2 = Mo/2 — yMi, which must 
be also hermit ian, given that both Mq and Mi have this 
property. 


Using that = —1, Eq. (|15|) can be rewritten as 


TyMiT = T, and this boundary condition becomes 


MT = T (17) 

with the definition M = TyMi. It follows from the sim¬ 
ilar block-structure of Ty and Mi, and from T’^'T = 1, 
that M is a diagonal matrix, the first d elements of which 
are 1 along the diagonal, and the last A — d are —1. From 
here, it is straightforward to see that {Jy,M} = 0, M 
is hermitian and unitary, and Mi must be unitary too. 
From {Ty,M} = 0, the properties of Ty and the basis- 
invariance of the trace, it also follows that M is traceless, 
therefore it has two eigenvalues of 1 and two of —1, mak¬ 
ing d = 2, and all the blocks in 2x2. Another 
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way to show this can be found in the basis where is 
a real antisymmetric matrix, and Mi is consequently a 
real symmetric matrix. The trace of the product of two 
such matrices is always 0. At this point, we see that M 
has two eigenvectors with eigenvalues of 1 , and two other 
eigenvectors with eigenvalues of —1. Note that the basis 
consisting of these eigenvectors is not fully specified as 
long as these pairs of eigenvectors can be rotated in their 
respective two-dimensional subspaces. We can specify 
the basis by fixing Y. For example, if T = cto, then 
Jry = —iT 2 ( 8 )cro and Mi = —ri ( 8 )<jo. Independently from 
the specific choice of T, we can call this basis the eigen- 
basis of M, because M = T 3 (g) do is diagonal. The sub¬ 
space of M with eigenvalues ±1 is mapped by into the 
subspace with eigenvalues which is a basis-invariant 
feature of course. 


From equations (13) and (16), we obtain that 
T’f'M 2 T = 0 , which means that the top left 2 x 2 block 
of the hermitian matrix M 2 must be 0 in the eigen- 
basis of M. In this basis, one can easily show that 
(M + 1 )M 2 (M + 1) =0, which also holds however in ev¬ 
ery other basis. From here, we obtain (M - 1 - 1 )M 2 T = 0, 
which combines with ( [l^ to the boundary condition for 
the derivative of the 4-spinor: 


(M + 1)AT' = 0. 


(18) 


This 8 X 8 Hamiltonian is commuting with the time re¬ 
versal operator 


T = -Sa, 0Ty0CFy C, (21) 


where C is the operator for complex conjugation. Upon 
time reversal, the real components of the k vector change 
sign, thus also changing the valley, and the electron spin 
flips. The edges on the other hand, do not couple dif¬ 
ferent electron spin states, so in the basis defined at 
(20), two blocks of M appear along the diagonal pos¬ 
sibly mixing states in different valleys but only with 
the same spin. Allowing that the parameter angles de¬ 
fined at are different for the two spin states, we re¬ 
quire that diag(M(/, (/:?, 0, 6 >), M(/', (/:?', 0', 6 >')) commutes 
with T, where diag denotes a (block-)diagonal matrix 
with the arguments along the diagonal. This can be 
solved with I' = —and ip' = ip^ (j)' = 0 , and 0' = 0^ 
and unless a ^ = a' ^ = 0, electron-hole symme¬ 

try or valley symmetry cannot be used to further reduce 
the number of the four free parameters of M, unlike in 
Ref. \28\. 


There is another symmetry of H that swaps the valleys 
and flips the spin, but does not change the sign of the 
momentum relative to the K points. This is represented 
by the operator 


Note that if A = 0, then ( |18[ ) is trivial, and we formally 
return to the solution for gapped graphene with the only 
boundary condition of Eq. 0 to be solved with four 
linearly independent functions at the edge. 

The two boundary conditions for T and T', Eqs. 0 
and (18) represent the first main results of this work. 


Both equations are expressed with the same M matrix, 
which is hermitian, unitary, and satisfies {Ty,M} = 0. 
These are the same conditions that led to Eq. (2.8) of 
Ref. ESI: 


M = sin / To (g) (^-cr) -j- cos / (iz-r) (g) (/la-ct), (19) 


where ly and fi are three-dimensional unit vectors, such 
that ^ and fj, are orthogonal to each other and also to rib, 
which is in turn normal to the boundary, and pointing 
outwards. The vector 1 / can be parametrized by the polar 
angle ip and the azimuthal angle 0, and the vectors jj, and 
^ with a single polar angle 0 and ^ -j- 7r/2. The optional 
minus sign in the latter, can be merged with the sign of 
the angle /. Assuming that no external magnetic field is 
present, we now invoke time reversal symmetry. In order 
to explicitly take into account the spin degree of freedom, 
we extend the Hilbert space. The Hamiltonian H must 
be extended with the same two blocks along the diagonal, 
but this time with 5 ^, a! ^ and P' parametrizing the 
K valley, and 5c, 5v, and (3 the K' valley. The 2x2 
blocks of H can now be labeled as 


- K'\ 

K't 

Ki 

K'i. 


( 20 ) 


Q = -5a; 0ry0CFyX (22) 

xdiag(—/c_, —/c_, —/c+, /c_, — k-)/k. 

This symmetry corresponds to a 7 r-rotation with respect 
to the y axis at an armchair edge (see Eig.[^ and results 
in spin-degeneracy there, in contrast to zigzag edges, 
where Q represents no symmetry and we observe spin¬ 
splitting. Unfortunately, this symmetry operation can¬ 
not be used to gain further information on M by requir¬ 
ing [M, Q] = 0, because of the momentum-operators in 
Q acting on edge states, which are superpositions of the 
bulk solutions with different k^. 

Apart from the mass term, and the valley-dependent 
quadratic terms a and /d, there are more differences be¬ 
tween graphene and the transition-metal dichalcogenide 
monolayers which we focus on here. While in graphene 
the basis functions are the orbitals in the two sublat¬ 
tices, the basis we are using here consists of the extended 
wave functions of the valence and conduction bands. This 
makes it less straightforward to see how the angles /, 99 , 
(j) and 0 relate to the geometry of a given edge. We will 
discuss here the two special cases of zigzag and armchair 
edges (Eig.j^an di, with a strong em phasi s on the zigzag 
edge for being energetically favorable^^^ and more rel¬ 
evant from the applications’ point of view. Erom Eig. 
it is clear that an armchair edge mixes the A" and K' 
valleys, while these remain non-equivalent in the ID Bril- 
louin zone of a zigzag edge. This suggests that the polar 
angle ip responsible for the valley mixing must be 0 or 
TT for zigzag, and from symmetry considerations, we also 
deduce that ip = tt{2 causing the maximal mixing of the 
valleys, must belong to the armchair orientation, just like 
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FIG. 3. Projections to the reciprocal spaces of zigzag and 
armchair nanoribbons. In armchair ribbons the valleys are 
fully mixed and all the K points are equivalent, while in zigzag 
ribbons, the K and K' points remain independent, the valley 
mixing is zero. 


in graphene.l^ This also means that (j) drops out of the 
equations (iT^ and (18) for zigzag, leaving only I and 0 to 


characterize the various zigzag edges. As it turns out, 0 
matters for armchair edges, but only if edge states from 
different edges can interfere. For a single armchair edge, 
it is also I and 0 that determine the edge states and their 
dispersion relation. Without studying the atomistic de¬ 
tails of an edge, based on symmetry considerations alone 
the values of I and 0 cannot be inferred within the contin¬ 
uum model. However, fitting to ab initio band structures 
is a straightforward way to relate the different edge ge¬ 
ometries to different values of these parameters. 

With ip = 0 for zigzag, M can be simplified to 


/ cos('i9+) sin('?^+) 0 0 \ 

sin('i9+) — cos('i9+) 0 0 I 

0 0 —cos{i^~) I ’ 

Vo 0 cos{h~) / 


(23) 

where = 0 ±1. The (p = tt case is already included if 
we shift '^9 by TT. Since the zigzag edge does not mix the 
valleys, and can be fitted separately. As we know 
from above, the same M^z also works for the opposite 
spin, but with —1. Note that I is expected to be small, 
for it is finite only because of the valley asymmetry. If 
all the parameters in H were the same in both valleys, 
then H would commute with —Ty 0 ayC (the 4x4 off- 
diagonal block of T), which in turn has to commute also 
with Mzz, resulting in / = 0 as in graphene (see Eq. (2.10) 
of Ref. [22 ). 


III. SINGLE EDGES AND NANORIBBONS 
A. General discussion 

Here we apply the M matrix method to various kinds 
of edges. Since we are interested in the in-gap edge states. 


the normal component of the wave vector k± will al¬ 
ways be imaginary. Out of the eight edge-solutions of the 
eigenvalue problem, we can only work with the four that 
are decaying away from the edge of a half-infinite system 
and we discard the diverging solutions. This leads to four 
unknown coefficients and eight equations defined by (17) 
and (18). However, the rank of each of the 4x4 matrices 
defining this linear system of equations is 2 , yielding in 
total four equations to solve for the four coefficients. The 
determinant of this 4x4 matrix is zero for the non-trivial 
solutions; this delivers the dispersion relation of the edge 
states, and the null-space of the matrix provides us with 
the coefficients. 

In a nanoribbon, the dispersion relation can be calcu¬ 
lated in the same way, but with eight states resulting in 
a determinant of an 8 x 8 matrix. Since the edge states 
are concentrated in the vicinity of the edges, in a few 
nm wide nanoribbon, they become independent. Further 
increments of the width do not affect the bands of the 
edge states any more, and the bands of the single edges 
independently add up to form the band structure of the 
nanoribbon. It is important to note that two edges, char¬ 
acterized by li ^2 and 0 i ^2 (or by i9f2 in zigzag), can be 
combined to form a nanoribbon, if we attribute —I and 
—0 (—'i^^) to the edge bordering the ribbon on the op¬ 
posite side relative to the one that was used to obtain 
these angles. Eg. in the case of an armchair ribbon with 
equivalent edges, /, 0, —I and —0 has to enter the above 
mentioned 8 x 8 matrix to give rise to the correct disper¬ 
sion relation originating from the combined effect of the 
two identical edges. 

This method can be generalized to ffakes with more 
complex shapes and mixed types of edges, like equilat¬ 
eral triangles, rectangles, hexagons, etc. Each edge has 
its respective 0 and /, and by studying the symmetries 
of the system of equations posed by the boundary condi¬ 
tions, it must be determined how to couple these to form 
the boundaries of the same flake. In the case of armchair 
edges, the relative distance and orientation of the edges 
also matter if the corresponding states are overlapping, 
meaning that the knowledge of the angles 0 is also re¬ 
quired for the accurate description. As we see at the end 
of this section through a simple example in graphene, this 
can be determined from an appropriate TB model. 


B. Application to M0S2 


To demonstrate the use of our method, we calculate the 
edge states in M 0 S 2 . To this end, we need the material 
parameters appearing in the Hamiltonian. Because of 
the mirror symmetry of the lattice to the x — y plane, H 
is block-diagonal in the ^-component of the real spin and 
takes the form of ( 20 ) with the valley-dependent sets of 


parameters {a, /3, 7 , 5v, and {a', /3', 7 , d? dl- We 
assume that no external held or special edge geometry 
breaks this symmetry of the bulk, and therefore the full 
Hamiltonian including M also preserves this symmetry. 
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Our continuum model, which describes the edge states 
as linear combinations of the bulk states using only the 
basis of the valence and the conduction bands, will not be 
able to describe all the in-gap states. For example, states 
that are almost entirely localized on the terminating row 
of sulfur atoms (see eg. Fig. 5 of Ref. [10]) cannot be 
expanded in the basis we use. The bulk wave functions 
of this basis require that the parameters appearing in H 
to be extracted from a bulk DFT study. The values of 
these bulk parameters are taken from Ref. [29], where 7 
is reported to be spin-independent: a = —2.16 eVA^, 
a' = -2.62 eVAA p = 4.09 eVA^ p' = 4.35 eVA^ 
7 = 2.76 eVA, Eg = e^-s^ = 1.821 eV, E'^ = = 

1.67 eV, Sc — = 3 meV and = 148 meV. Note 

that the contribution of the free Hamiltonian / (2me) = 
3.81 eVA^ is included in the a, a', [3 and P' parameters, 
where tUq is the free electron mass. 

In what follows, we fit and compare our results to the 
DFTB study of Erdogan et al for nanoribbonsAl¬ 
though SOI was not considered, the dispersion curves of 
sufficiently wide (32.8 A) nanoribbons with several dif¬ 
ferent kinds of edges were calculated, which offers the 
opportunity to fit'd edge-by-edge assuming that the two 
sides are independent. Such DFT calculations are typ¬ 
ically done with nanoribbons and not with single edges 
because of the periodic boundary condition. However, we 
can consider the edge bands independent in nanoribbons 
over a few nm wide, because exponentially decaying edge 
states on either sides become negligible at the opposite 
edge and they do not interfere. Here we first adopt the 
assumption that the DFTB band structure is actually 
spin-resolved and belongs to the spin-up states in the 
K' valley, and then we repeat the calculation with the 
assumption that the same bands belong to spin-down 
states in the same valley. With this fitting to spin- 
independent bands, we aim to demonstrate the use of 
our method, rather than providing quantitatively accu¬ 
rate results. Accuracy could be obtained from fitting 
to non-magnetic, spin-resolved bands, but to our best 
knowledge, no such study has been published to this day. 

The fitted spin-up bands are plotted in Fig.|^ together 
with the original DFTB bands and nanoribbon struc¬ 
tures from Ref. EBP To fix the energy offset of the 
edge bands, we match the center of the bulk gap to the 
center of the nanoribbon energy gap by a rigid energy 
shift for the whole band structure. Then we adjust the 
energy of our edge band at the K' point to that of the 
DFT-calculated band by tuning the only free parameter 
'd for each single edge band separately. The fitted'd val¬ 
ues are unique apart from the 27r periodicity, and these 
were used for the plots of Fig. i as well.l^ It is apparent 
from Fig. that all four band structures share an in-gap 
state with a negative slope. This suggests that this band 
belongs to the same edge that all the four nanoribbons 
have in common, namely the purely sulfur terminated 
zigzag edge (the top edges of the ribbons in Fig. [^. To 
this band, we can fit the angle with almost the same 
value in every cases. The edges on the other side are 


described by '^2 , and in this case, we can only fit to 
one band without ambiguity even though in two cases 
(Fig. B a) and (c)), other in-gap bands appear in DFT, 
which are missing from the continuum model. Note that 
the edge with one S atom per Mo supports no edge state 
at all, see Fig. 


The deviations from the DFTB band structure can be 
explained with the following factors: (i) in the DFTB 
calculation, the spin-orbit coupling was neglected; (ii) 
the input parameters for our calculation (a, and 

Eg) were extracted from a bulk DFT calculatiorP^ with 
different exchange-correlation potential, etc; and finally 
(hi) we used an effective, 2-band Hamiltonian neglecting 
the trigonal warping and the cubic terms.Regarding all 
these differences, the agreement between the band struc¬ 
tures is reasonably good. The most significant deviation 
manifests itself in the larger steepness of the bands that 
have lower energy at the K' point. Note however that 
this band is moderately steeper also in Refs. [34] relative 
to our reference,!^ and the derivative of our bands agrees 
especially well with the tight binding results in Fig. 9 
from Ref. [26], even if we need to use slightly different'd 
angles to fit our bands to theirs. 

The corresponding components of the wave function 
can be seen in Fig.[^ where the real-valued coefficients of 
and p)c are plotted along the dimension normal to the 
various kinds of zigzag edges in Fig.[^ Here, the left edge 
always has the same, purely sulfur terminated structure, 
whereas the right edge has different reconstruction in the 
four cases. Since the S-half ribbon has no edge states on 
the right (see Fig. 4d), the dotted lines only appear on 
the left. The differences in the wave function on the right 
manifest in slightly different values for the otherwise 
identical left edges (see Fig. [^, and here, this prevents 
the ipy curves and also the curves from exact overlap. 
This is due to the small, yet finite contribution of the 
electrons from the opposite edge 32.8 A away. On the 
left, the overall predominance of the conduction band 
can be seen, but on the right, the contribution of the two 
bands is roughly the same. Accordingly, the left edge 
state will always have the higher energy, see Fig.[^ Using 
Eq. (13), we can verify a direct relation between these 
curves, namely that -h 'ipcP'i^c — ^ holds at both 

edges. The ratio of the a' and P' parameters explains 
the relative height and steepness of the and p)c curves 
at the two ends. 


Note that in graphene, there are basically only two 
kinds of zigzag edges, and even they are related by sym¬ 
metry. One edge is always terminated by atoms of sub¬ 
lattice A, and the other one by atoms of sublattice B. 
In our case, these atoms are different elements, we lose 
inversion symmetry, and also the fact that molybdenum 
atoms form bonds with six out-of-plane sulfur atoms in¬ 
creases the number of possible edge terminations. In the 
case of the antisymmetric edges (Eig. [^, the purely S 
and purely Mo terminated edges correspond to the solid 
curves in Eig. on the left and on the right respectively. 
This S-Mo antisymmetry constitutes the major differ- 
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(a) Antisymmetric edges 



(b) S 2 -stripe edges 



(c) S-dimer edges 



(d) S-half edges 



= 1.795, =77 + 1.762. 



= 1.768, =77 + 1.840. 



= 1.791, =77 + 1.810. 



k [7/4] 

= 1.780. 


FIG. 4. On the left, various 12 dimers (32.8 A) wide zigzag nanoribbons from Ref. [21], with the corresponding band structures 
on the right. The edge state bands from the continuum model for spin-up states are denoted by thick, green (gray) lines, and 
those from the DFTB calculation from Ref. [21] are denoted by thin, black lines. The dotted lines show the K' point, and the 
wave number is measured from the F point in these figures. The thick green (gray) curves correspond to the states at the two 
edges fitted with 'd^ and 7^^, where = +2 = 0. The dots represent the data points the curves were fitted to. 





















































FIG. 5. Coefficients of 'ijj-y and the wave function at 

kx = 0 in the Kl and K'^ valleys at the zigzag edges, within 
the width of two dimers y/Sa. Solid lines describe the an- 
tisy mme tric edges (Fig. |4a| ), dashed lines the S 2 -stripe edges 
(Fig.|4b|) and dashed-dotted lines the S-dimer edges (Fig. [4^. 
Dotted lines appearing only on the left describe the S-hal f rib - 
bon, which does not support edge states on the right (Fig. 4d). 
The wave function is calculated for independent edges with 
the'd values from Fig. Note that all the components of T 
are real, and the maximum of '0c (^v) on the left (right) is 
always normalized to 1. For the corresponding energies, see 
the band structures in Fig. at the K' point. 


ence between the wave functions on the two sides. The 
dangling bonds of the Mo atoms offer further possibility 
for various reconstructions of the edge with additional S 
atoms thus modifying the picture to some lesser extent on 
the right, but it does not change qualitatively. As we see, 
the obtained values for the angles on the right {^ 2 ) 
and on the left are also forming two distinct groups, 
emphasizing that the sublattice asymmetry is more im¬ 
portant than the effect of the various edge reconstruc¬ 
tions on the Mo-terminated side. This is also underlined 
regarding the energy of these states: the maximal energy 
difference nt kx = 0 between the two sides is 531 meV 
and the slopes of the bands differ in sign. The maximal 
energy difference between the various reconstructions on 
the right relative to the middle of the gap is 351 meV, 
and the slopes of the bands are almost the same. And 
finally, the maximal difference on the left due to edge 
states on the right reaching the left edge is only 26 meV. 
Remember that in the context of the continuum model, 
the different values of can be understood to be due to 
the small, yet finite contribution of the states from the 
opposite edges at a distance of 32.8 A. Otherwise, one 
kind of edge is characterized by one pair of angles 
which in the case of the purely S terminated left edge 


would be around 1.78. The precise determination of 
would necessitate fitting to ab initio calculations with 
somewhat wider ribbons, or the simultaneous fitting of 
i9f and 1 ^^ in a ribbon explicitly considering the interfer¬ 
ence of the edge states. This accuracy however would be 
justified only if we were fitting to a spin-resolved band 
structure accounting for the spin-orbit interaction as well. 

There is one more detail to zigzag M 0 S 2 nanoribbons 
which may be potentially important. Namely, experi¬ 
ments have shown e vidence of weak ferromagnetism in 
M 0 S 2 nanostructures.^ ^ * ^^ * ^^ * It was suggested that this 
can be explained by recent calculations on zigzag M 0 S 2 
monolayer nanoribbons, which indicate that the ground 
state of these nanoribbo ns can be ferromagnetic, depend¬ 
ing on the tensile strain,I^ ^J^ and also on the edge passi¬ 
vation (eg. H-saturation) The consequent energy 

splitting of the edge states, the magnetization splitting, 
depends on the average magnetic field felt by the electron 
in a given state. This is mainly determined by relation 
of the decay length of the state to the magnetic profile of 
the edge. Around the K points, it is difficult to separate 
the effect of the SOI and that of the edge magnetization 
(EM). However, by studying the edge-b ands in the rel¬ 
evant literature on DFT band structureswe can 
see that the combined effect of the two is more or less a 
rigid shift along the energy axis for the two edge-bands 
that can be described by our model. This shift depends 
on the band, but there seems to be negligible dependence 
on k within one band, and so the change in the associ¬ 
ated group velocities at the K points are also negligibly 
small. Therefore, we assume that around the K points, 
where the effective Hamiltonian H is accurate, the effect 
of EM on the one-dimensional edge-states can also be ac¬ 
counted for by a simple energy shift of the bands. To our 
best knowledge, the band structure of the non-magnetic 
state has not been published yet for such edges, although 
calculations were made eg. for the total energy difference 
of the magnetic and non-magnetic states (see Eig. 6 of 
Ref. [24]). On average, this reaches only tens of meV, 
but some bands are affected by the EM more than oth¬ 
ers. This input would be needed to perform a proper 
fitting for i^^ 2 ^ after which, a simple energy-shift of the 
bands could restore the DET band structures in a good 
approximation. Eig. shows the bands at the purely 
sulfur terminated zigzag edge in the non-magnetic state, 
with the assumption that = 1.795. The split¬ 

ting at /c = 0 is 38 meV. The order of magnitude of this 
energy gap can be estimated from the spin-splitting of 
the states '0v and '0c, and the ratio of their contribution 
to the given edge state (left part of Eig. |^. This splitting 
is significantly smaller than t he values over 300 meV re¬ 
ported in the literaturej ^^ l ^^ l This could suggest that the 
EM is more important than the SOI for these states, 
provided that the above equality of the angles holds. 
We find however that the gap changes by 1 meV for a 
change of 0.001 in the value of I around 0 ^ 1.795, and 
for example, / = 0.2 would already result in a splitting 
of 343 meV due only to SOI in good agreement with the 
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FIG. 6. Band structure of the purely sulfur terminated zigzag 
edge in M 0 S 2 (top edges in Fig.[^, where electrons are prop¬ 
agating in opposite directions in the K and K' valleys. The 
bands are calculated for the non-magnetic state, that is the 
spin-splitting is only due to the bulk SOI, where the assump¬ 
tion = 1.795 is used. Down-arrows and dashed 

lines denote down-spin states, up-arrows and solid lines up- 
spin states. The left-propagating states in the K valley are in 
black, and the right-propagating states in the K' valley are 
in green (gray). 


DFT calculations. Nevertheless, to know exactly the rel¬ 
ative strength of these two effects, together with the pre¬ 
cise values of a fit to a non-magnetic band structure 
would be needed. Note that this large spin-splitting com¬ 
bined with the fact that K- and iF'-electrons are prop¬ 
agating in opposite directions could be used to generate 
spin- and valley-polarized edge currents at room tem¬ 
perature and above .1^ The EM can further enhance the 
spin-splitting for one of the valleys. In a first approxi¬ 
mation, the magnetization is expected to shift the 
band up (or down) by the same amount as the K'l band 
shifts down (or up). A similar statement holds for the 
other two bands as well, but the magnetization splitting 
of these will be slightly different due to the valley asym¬ 
metry. 

In Ref. [ 34 ] , the edge states of M 0 S 2 nanoribbons are 
studied in relation to the width-dependent electrostatic 
polarization, with the help of first-principles DFT calcu¬ 
lations. According to this study, the wider the ribbon 
is, the more these polarization charges are screened, but 
in the case of incomplete screening in narrower ribbons, 
a transverse electric field is induced by the polarization 
charges. As a result, a potential energy difference ap¬ 
pears between the edges, and the edge bands are pushed 
upward at one side of the ribbon, and pushed downward 
at the other side. This phenomenon is independent of 
the spin state, but otherwise, in our model, it can be ac¬ 
counted for in a similar fashion as EM, with a rigid shift 
of the bands. _ 

Recently, several workP^^ used a “gapped-graphene” 
type k p Hamiltonian to model the bulk electronic struc¬ 
ture of monolayer M 0 S 2 in the vicinity of the K {K') 
points. While such analogy to graphene can be use¬ 


ful to understand bulk properties, the analogy breaks 
down at the edges. To highlight the difference to gapped 
graphene, it is instructive to consider a zigzag edge in the 
limit of vanishing or negligible a and (3 parameters and 
equal energy gaps for the two valleys and spins. Here 
we assume that the restored electron-hole symmetry of 
the bulk is not broken by the edge. Erom this, it fol¬ 
lows that 9 = = 7TI2^ and we obtain edge states with 

Dirac-like, linear dispersion 5 = for bands belong¬ 

ing to the K and K' valleys respectively. The normal 
component k± = iE^/ ( 27 ) depends only on material pa¬ 
rameters. These results are independent from the edge 
geometry, that is all zigzag edges are equivalent in the 
model as long as a = = 0. Differences between the 

edges can be taken into account only if finite values for a 
and (3 are allowed. But why can we not find these in-gap 
states in gapped graphene^^ To answer this, we recall 
that the basis of graphene is a sublattice basis, in which 
the edge inherently breaks the electron-hole symmetry, 
since all the terminating sites belong to one particular 
sublattice, and all have the same on-site energy, which is 
also the energy of the fiat bands (see Eig. 1 of Ref. m)- 
In M 0 S 2 , the basis functions are extended wave functions 
of the valence- and conduction-bands, and the terminat¬ 
ing sites do not belong to either of these. Therefore, it is 
conceivable that none of the two basis functions is given 
preference by the edge, thus not causing any symmetry 
breaking. Based on these findings, we expect the pres¬ 
ence of edge states showing the preserved electron-hole 
symmetry of the gapped system. 

In the case of an armchair ribbon, provided there is no 
asymmetric edge reconstruction, we have the same kind 
of edges on both sides: li = I 2 and 0i = ^ 2 - We have 
two different linear combinations of states localized to the 
opposite edges, and these combinations are split in en¬ 
ergy, but they are degenerate in spin. In narrow ribbons, 
where the edge states do not vanish completely at the op¬ 
posite edges, the difference Acj) = 02 — ^1 affects the band 
structure, but the angles 01^2 individually do not. In wide 
ribbons, where the opposite edges do not interfere, the 
angles 0 have no effect on the dispersion, the splitting 
vanishes, and we end up with fourfold degenerate edge 
states: twice for the spin and twice for the two linear com¬ 
binations, where the latter ones can eventually be split 
into states belonging to their respective edges only. Here, 
the angles 0 turn out to be arbitrary phase differences be¬ 
tween the valleys, as in ~i^c 6 xp(i 0 ), exp(i 0 )}, 

and this phase difference does not change upon time re¬ 
versal. 

It is apparent from Eig. 3(b) of Ref. [ 21 ] that the band 
structure of the edge states does not change qualitatively 
with the width in the given range, and the valley-mixing 
effect of A0 can be neglected. Because of the size of our 
basis, we cannot describe all the bands, but only one of 
the degenerate flat bands around the iF* point (Eig. [^, 
denoted by F in Fig. 3(b) of Ref. [21]. For small values of 
/, almost in the full range of we see a flat band indeed 
that could be matched with any of the two between 0 
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and 1 eV. Without further input however, eg. regarding 
the composition of these states, we cannot tell which one 
to fit to. 

Note however that in graphene nanoribbons, we can 
also describe extended edge-to-edge bulk states with the 
M matrix method, where Acj) also has an effect on the 
dispersion of these states. This can be used to determine 
its value eg. by means of a TB model, which enables 
the analytical determination of band structure, like it 
has been done in Ref. For the spectrum at the K 


point, these authors found Sr = t 


1 + 2 cos 


A^d + l 


, where 


t is the transfer integral between nearest-neighbor car¬ 
bon sites, and r G { 1 ... A^d} is an integer indexing the 


A^d + l 


27r 

3 


bands. For the lowest, positive energy Sr 
is minimal with r = tq. It can also be shown easily 
that Acj) = ±2Werol"f holds, where W = and 

7 = ^ta. Combining all this, we obtain that Ac/) = ±|7r 
if A^d modulo 3 is 0 or 1, and A0 = 0 if A^d modulo 3 
equals 2 . The sign of Acj) is irrelevant for armchair 
nanoribbons because of their symmetric edges, this is 
why we cannot find a definite sign from the energy spac¬ 
ing. Nevertheless, it is sensible to think of this effect 
as if every additional dimer added 27r/3 to 02 counting 
from the narrowest ribbon with A^d = 3, 0i = 0 and 
02 = 27r/3. As for the in-gap states in TMDC armchair 
nanoribbons, we speculate that A0 may be determined 
by the width of the ribbon as well. The main arguments 
are that the individual angles 0i and 02 do not affect the 
in-gap bands, but only their difference, and the effect of 
A0 vanishes with the vanishing overlap of opposite lying 
edge states in an infinitely wide ribbon. These findings 
suggest that the relative geometry of the edges matter, 
as in the case of graphene, where A0 is directly pro¬ 
portional to the ribbon width W. The determination of 
the exact relation between A0 and W in TMDCs is an 
open probl em, w hich could probably be solved by a TB 
calculation 


monolayers within a continuum model. We used an ef¬ 
fective, two-band, k p Hamiltonian around the K point, 
where the electron-hole symmetry breaking quadratic 
terms were taken into account. We modeled the effect 
of the edges with the M matrix method, and after re¬ 
ducing the number of free parameters by symmetry con¬ 
siderations, we found that different edge geometries can 
be described with two or three scalar parameters varying 
between 0 and 27r. We discussed the case of single edges 
and interacting edges of nanoribbons, and we fitted the 
free parameters to ab initio band structures. Focusing 
mainly on zigzag edges, we analyzed the edge states and 
their dispersion relation in M 0 S 2 in particular, and we 
demonstrated a good agreement with DFT based tight 
binding calculations. 

Note that our method is not restricted to in-gap states, 
but it can be applied to low-energy bulk states as well, 
as long as they are close to the conduction band or 
the valence band edge. Therefore our method can be 
used to study bulk band states in edge-terminated nano¬ 
structures, such as quantum dots and nanoribbons 
This also opens up the possibility to fit the angles 
for the two edges using these extended bulk states. The 
advantage of this would be that the bulk bands are not 
much affected by the EM in wider ribbons, but the dis¬ 
advantage is that the bulk states connect the edges and 
a simultaneous fitting for i9f and is needed. The 
ribbons cannot be too wide either so that the quantized 
bulk bands are well-separated in energy. Alternatively, 
as demonstrated in our work, fitting to non-interacting 
in-gap edge states is the simplest way to obtain i9f2i but 
for accuracy, non-magnetic band structures need to be 
calculated first. 
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